clear;

V=-1;
doping=6.25;

figure('position',[300,300,200,500]); hold on; box on;
ii=0;

for EE=-0.9:0.1:-0.1
t=0.67;
U=8;
rootPath = ['extHubbard_U8_V',num2str(V),'_doping',num2str(doping),'/'];
cptSpectrum = load([rootPath,'AkwCPTSpectra.txt']);
enrgList = load([rootPath,'enrgList.txt']);
momList = 1:size(cptSpectrum,1);
if size(cptSpectrum,1) == 501
    momList = [1:200 200+sqrt(2)*(1:100) (200+sqrt(2)*100)+(1:201)];
end
specInfo = load([rootPath,'specInfo.txt']); mu = specInfo(4);

cptSpectrumG=cptSpectrum;

%Fermi Dirac function
kb=8.6e-5;
FD=1./(exp((enrgList-mu)*t/(kb*20))+1);
GFD=normpdf(enrgList*t,median(enrgList*t),0.1);
FDG=median(diff(enrgList*t))*conv(GFD,FD,'same');


G=normpdf(enrgList*t,median(enrgList*t),0.18);
for i=1:length(momList)
    cptSpectrumG(i,:)=median(diff(enrgList*t))*conv(G,cptSpectrumG(i,:),'same').*FDG; 
end
SpectrumGSym=[fliplr(cptSpectrumG'),cptSpectrumG(2:end,:)'];
S=SpectrumGSym;

E=(enrgList-mu)*t;

for iii=1:length(E)
    if E(iii)>EE
        break;
    end
end
    
mdc=S(iii,:);
mdc=mdc./max(mdc)*0.5;

k=linspace(1,length(mdc),length(mdc));
k=k-mean(k);
k=k/max(k)*0.788;

offset=3.3/3-1.1;%y axis offset
inteval=0.53/2;%y axis offset
plot(k,mdc+inteval*ii+offset,'r');


ii=ii+1;

end


xlim([0.1,0.788]);
ylim([0,2.7])
xlabel('k_{||} (�)');
ylabel('intensity (arb. units)');
set(gca,'PlotBoxAspectRatio',[0.4 3 1],'FontSize',12);
